Projekt z zakresu analizy danych - wpływ rozkładu danych na MSE regresji

Author

Diana Morzhak, Mateusz Walo, Dominika Zydorczyk

Wstęp i krótkie streszczenie projektu

Celem projektu jest zbadanie, w jaki sposób rozkład zmiennych objaśniających wpływa na wartość błędu średniokwadratowego ( \(MSE\) ) w modelu regresji liniowej.

Zagadnienie to jest istotne zarówno z teoretycznego, jak i praktycznego punktu widzenia - w szczególności w uczeniu maszynowym i analizie danych, gdzie poprawność założeń dotyczących rozkładu predyktorów wpływa na stabilność i trafność estymacji parametrów modelu.

Rozpatrujemy wpływ rozkładu zmiennych objaśniającyh na błąd średniokwadratowy (MSE) w następującym modelu regresji liniowej:

\[Y = 0.5 \ + \ X_1 \ + \ 0.6X_2^2 \ + \ \epsilon \]

gdzie \(X_1\) i \(X_2\) są zmiennymi losowymi pochodzącymi z rozkładów zależnych od rozpatrywanej sytuacji (patrz niżej), natomiast \(\epsilon\) jest szumem pochodzącym z rozkładu \(N(0,0.1)\) - jest tak, ponieważ chcemy się skupić wyłącznie na wpływie rozkładu predyktorów. MSE definiujemy następująco:

\[ MSE=\frac{1}{n} \sum_{i=1}^n(Y_i-\hat{Y_i})^2 \]

gdzie \(Y=(Y_1, Y_2, ..., Y_n)\) jest wektorem wartości obserwowanych, natomiast \(\hat{Y}=(Y_1, Y_2, ..., Y_n)\) jest wektorem wartości przewidywanych przez model, a \(n\) oznacza liczbę obserwacji w zbiorze danych. Dla całego projektu przyjmujemy, że \(n = 120\) , co pozwoli nam na bezpośrednie porównywanie uzyskiwanych wyników.

Lista rozpatrywanych przypadków

Poniżej przedstawiono wszystkie przypadki, w których generowane są zmienne losowe \(X_1\) i \(X_2\) . Opisano zarówno przypadki niezależne, jak i zależne, z uwzględnieniem możliwej zależności liniowej i nieliniowej tych zmiennych oraz wpływu zmiennej ukrytej na nie:

  1. double_snorm
    W tym przypadku zarówno \(X_1\) , jak i \(X_2\) pochodzą z rozkładu normalnego standaryzowanego: \(X_1\sim \mathcal{N}(0,1), \ X_2 \sim \mathcal{N}(0,1)\)

    Ten przypadek jest naszym punktem odniesienia, ponieważ normalny rozkład standaryzowany jest symetryczny oraz ma umiarkowany rozrzut danych.

  2. double_unif_narrow
    Zarówno \(X_1\) , jak i \(X_2\) są losowane z wąskiego rozkładu jednostajnego, gdzie zakres danych jest ograniczony, a gęstości prawdopodobieństwa sa równomierne: \(X_1 \sim \mathrm{U}(-3,3), X_2 \sim \mathrm{U}(-3,3)\)

  3. double_unif_wide
    W tym przypadku zmienne są losowane z szerszego rozkładu jednostajnego: \(X_1\sim \mathrm{U}(-5,5), X_2 \sim \mathrm{U}(-5,5)\)

  4. norm_exp
    Zmienna \(X_1\) ma rozkład normalny z większym odchyleniem standardowym, a \(X_2\) pochodzi z rozkładu wykładniczego, co poszerza nasze badania o rozkład asymetryczny: \(X_1 \sim \mathcal{N}(0, 2.9), \quad X_2 \sim \mathrm{Exp}(\lambda = 1.2)\)

  5. snorm_exp
    Tutaj \(X_1\) ma rozkład normalny standaryzowany, a \(X_2\) - wykładniczy. Przypadek ten w zestawieniu z przypadkiem 4. pozwala zobaczyć wpływ wariancji rozkładu normalnego na MSE: \(X_1 \sim \mathcal{N}(0,1), \quad X_2 \sim \mathrm{Exp}(\lambda = 1.2)\)

  6. bimodal_snorm
    W tym przypadku \(X_1\) ma rozkład bimodalny, a \(X_2\) - normalny standaryzowany: \(X_1 \sim 0.43\cdot \mathcal{N}(-3,0.8) + 0.57\cdot \mathcal{N}(3,0.8), \quad X_2 \sim \mathcal{N}(0,1)\)

    Dzięki temu zobaczymy, jak zmienia się MSE w przypadku gdy mamy dwie mody w jednym z rozkładów.

  7. exp_norm
    Zmienna \(X_1\) pochodzi z rozkładu wykładniczego, a \(X_2\) z normalnego: \(X_1 \sim \mathrm{Exp}(\lambda = 1.2), \quad X_2 \sim \mathcal{N}(0, 2.9)\)

    Tutaj sprawdzamy sytuację w pewnym sensie odwrotną do przypadku 4., czyli co się dzieje, gdy to zmienna \(X_1\) ma duży ogon rozkładu.

  8. t_snorm
    Zmienna \(X_1\) pochodzi z rozkładu t-Studenta z 3 stopniami swobody, a \(X_2\) z rozkładu normalnego standaryzowanego: \(X_1 \sim t(df=3), \quad X_2 \sim \mathcal{N}(0,1)\)

  9. t_norm
    Zmienna \(X_1\) pochodzi z rozkładu t-Studenta, a \(X_2\) z normalnego rozkładu o większej wariancji: \(X_1 \sim t(df=3), \quad X_2 \sim \mathcal{N}(0,2.9)\)

  10. t_exp
    Zmienna \(X_1\) pochodzi z rozkładu t-Studenta, natomiast \(X_2\) z rozkładu wykładniczego: \(X_1 \sim t(df=3), \quad X_2 \sim \mathrm{Exp}(\lambda = 1.2)\)

W przypadkach 10.-12. analizujemy wpływ grubych ogonów i wartości odstających na MSE. Te trzy sytuacje możemy porównywać z dowolnym z przypadków 1-9, gdzie zmienna objaśniająca \(X_2\) ma dokładnie taki sam rozkład jak w przypadkach 10-12.

  1. Korelacja liniowa
    linear_corr_snorm Obie zmienne pochodzą z rozkładu normalnego standaryzowanego, ale wprowadzona jest między nimi korelacja \(\rho\) : \((X_1, X_2) \sim \mathcal{N}\Big(\mu =\begin{pmatrix}0\\0\end{pmatrix},\Sigma =\begin{pmatrix}1 & \rho \\ \rho & 1\end{pmatrix}\Big)\)
    Taki przypadek pozwala badać wpływ liniowej zależności między predyktorami

  2. Zależność nieliniowa nonlinear_corr Zmienna \(X_2\) pochodzi z rozkładu normalnego standaryzowanego. Zmienna \(X_1\) jest nieliniową funkcją zmiennej \(X_2\) z dodatkiem szumu \(\eta\) : \(X_2 \sim \mathcal{N}(0,1), \quad X_1 = a \cdot e^{b X_2} + \eta, \quad \eta \sim \mathcal{N}(0,0.1)\)
    Pozwala to na modelowanie przypadków, w których zależność między predyktorami jest wyraźnie nieliniowa

  3. Zmienne zależne od zmiennej ukrytej Z hidden_impact

Tutaj zarówno \(X_1\) , jak i \(X_2\) zależą od wspólnej zmiennej ukrytej \(Z\) : \(X_1 = 2.7 Z + 12, \quad X_2 = 1.4 Z + 15, \quad Z \sim \mathcal{N}(0,1)\)
Taki przypadek umożliwia zbadanie wpływu zmiennej ukrytej na obserwowane zmienne objaśniające

Razem zbadamy 17 przypadków. Dla każdego z nich generować będziemy po 200 ramek danych, co da nam wystarczającą próbkę do porównywania wyników. Dla każdej z utworzonych ramek danych będziemy dopasowywać model regresji liniowej, zapiszemy uzyskane wyniki i porównamy je. Tym samym wyciągniemy najciekawsze informacje z wpływu różnych rozkładów predyktorów na wartość błędu średniokwadratowego tworzonej regresji.

Generowanie danych

Funkcja jako narzędzie

Wprowadzamy funkcję Y_true , która oblicza wartości wektora Y na podstawie wygenerowanych danych. Ponadto dodajemy epsilon oraz eta, czyli szum, aby nadać prawdziwym wartościom \(Y\) losowości. W ten sposób tworzymy “prawdziwe, rzeczywiste” dane.

Do generowania danych do każdego z wymienionych we wstępie przypadków rozkładu danych zmiennych objaśniających wykorzystujemy funkcję generate_data . Uwzględnia ona wszystkie założenia rozpatrywanych przez nas w tym projekcie scenariuszy. Generuje odpowiednie dane w zależności od wprowadzonego argumentu case (co tłumaczymy jako ‘przypadek’). Następnie wygenerowane dane są zapisywane do ramki danych, która zawiera 4 kolumny:

  • X1 - wartości pierwszej zmiennej objaśniającej

  • X2 - wartości drugiej zmiennej objaśniającej

  • Y - obliczone na podstawie wzoru (numer wzoru) wartości Y

  • Eps - biały szum

Funkcja zwraca tę ramkę danych oraz meta dane generatora rozkładu, czyli listę parametrów użytych podczas generowania danych. Parametry te można poznać poniżej. A zatem wynik funkcji generate_data to lista posiadająca dwa elementy:

  1. df - wygenerowana ramka danych zgodnie z wybranym przypadkiem case
  2. meta - meta dane przechowujące wartości parametrów użytych do generowania danych
Pokaż/Ukryj kod 
# Funkcja służąca do obliczenia prawdziwych wartości Y
Y_true <- function(x1, x2, epsilon){
  return (0.5+x1+0.6*x2*x2+epsilon)
}
Pokaż/Ukryj kod 
generate_data <- function(case, 
                          n = 120, # liczba generowanych obserwacji
                          rho = NULL, # współczynnik korelacji (przypadek zależności liniowej)
                          lambda = 2, # lambda dla rozkładu Poisson'a
                          exp_rate = 1.2, # intensywność dla rozkładu wykładniczego
                          a = NULL, # wartość pierwszego współczynnika dla zależności nieliniowej
                          b = NULL # wartość drugiego współczynnika dla zależności nieliniowej
                          ){
  # definicja szumu
  epsilon <- rnorm(120, sd = 0.1)
  # Globalna definicja szumu dla przypadku zmiennej zależnej
  eta <- rnorm(120, sd = 0.1)
  
# sprawdzamy, czy został podany obowiązkowy parametr dla danego przypadku generowania danych
  needed <- function(param, name){
    if(is.null(param)){
      stop(sprintf("Przypadek '%s' wymaga argumentu '%s'", case, name))
    }
  }
  
  # przypadek, gdy oba predyktory mają rozkład normalny standaryzowany
  if(case=="double_snorm"){
    x1 <- rnorm(n)
    x2 <- rnorm(n)
  # Zmienne X_1 i X_2 pochodzą z rozkładu jednostajnego o parametrach (-3,3)
  }else if(case == "double_unif_narrow"){
    x1 <- runif(n, min = -3, max = 3)
    x2 <- runif(n, min = -3, max = 3)
  # Zmienne X_1 i X_2 pochodzą z rozkładu jednostajnego o parametrach (-5,5)
  }else if(case == "double_unif_wide"){
    x1 <- runif(n, min = -5, max = 5)
    x2 <- runif(n, min = -5, max = 5)
  # Zmienna X_1 pochodzi z rozkładu normalnego z większym odchyleniem standardowym,
  # natomiast zmienna X_2 pochodzi z rozkładu wykładniczego z lambdą = 1.2 (parametr intensywności)
  }else if(case == "norm_exp"){
    needed(exp_rate,"exp_rate")
    x1 <- rnorm(n, sd = 2.9)
    x2 <- rexp(n,rate = exp_rate)
  # Zmienna X_1 pochodzi z rozkładu normalnego standaryzowanego, a zmienna X_2
  # z rozkładu wykładniczego o lambdzie równej 1.2
  }else if(case == "snorm_exp"){
    needed(exp_rate,"exp_rate")
    x1 <- rnorm(n)
    x2 <- rexp(n, rate = exp_rate)
  # Zmienna X_1 pochodzi z rozkładu bimodalnego, a zmienna X_2 z normalnego
  # standaryzowanego
  }else if(case == "bimodal_snorm"){
    bimod_mean1 <- -3
    bimod_mean2 <- 3
    bimod_sd <- 0.8
    p <- 0.43
    u <- runif(n)
    x1 <- ifelse(u<p, rnorm(n, bimod_mean1, bimod_sd), 
                      rnorm(n,bimod_mean2,bimod_sd))
    x2 <- rnorm(n)
  #Przypadek, gdy zmienna X_1 pochodzi z rozkładu wykładniczego,
  # a zmienna X_2 z rozkładu normalnego
  }else if(case == "exp_norm"){
    needed(exp_rate,"exp_rate")
    x1 <- rexp(n, rate = exp_rate)
    x2 <- rnorm(n, sd = 2.9)
  }
  #X_1 pochodzi z rozkładu t-Studenta o trzech stopniach swobody,
  # a X_2 z rozkładu normalnego standaryzowanego
  else if(case == "t_snorm"){
    x1 <- rt(n,df = 3)
    x2 <- rnorm(n)
  #X_1 pochodzi z rozkładu t-Studenta o trzech stopniach swobody,
  # a X_2 z rozkładu normalnego
  }else if(case == "t_norm"){
    x1 <- rt(n, df = 3)
    x2 <- rnorm(n,sd = 2.9)
  #X_1 pochodzi z rozkładu t-Studenta o trzech stopniach swobody,
  # a X_2 z rozkładu wykładniczego
  }else if(case == "t_exp"){
    needed(exp_rate,"exp_rate")
    x1 <- rt(n, df = 3)
    x2 <- rexp(n, rate = exp_rate)
  # Zmienne losowe X_1 i X_2 są liniowo skorelowane i pochodzą obie 
  # z rozkładu normalnego standaryzowanego
  }else if(case == "linear_corr_snorm"){
    needed(rho,"rho")
    sigma <- rbind(c(1,rho), c(rho,1)) #macierz kowariancji
    mu <- c(0,0) #średnie rozkładów
  # Zmienna X_1 jest nieliniowo zależna od zmiennej X_2
  # a zmienna X_2 pochodzi z rozkładu normalnego standaryzowanego
  }else if(case == "nonlinear_corr"){
    needed(a,"a")
    needed(b,"b")
    x2 <- rnorm(n)
    x1 <- a*exp(b*x2)+eta
  }else if(case == "hidden_impact"){
    z <- rnorm(n)
    x1 <- 2.7*z + 12
    x2 <- 1.4*z + 15
  }
  
  meta <- list(case = case, n = n, 
               params = list(rho = rho, lambda = lambda, exp_rate = exp_rate))
  
  
  if(case == "linear_corr_snorm"){
    df <- as.data.frame(mvrnorm(n=n, mu=mu, Sigma=sigma))
    df <- df %>% 
    mutate(Y=Y_true(V1, V2, epsilon),
           Eps = epsilon)
    colnames(df) <-  c("X1","X2","Y", "Eps")
    return(list(df=df,meta=meta))
  }else if(case == "hidden_impact"){
    df <- data.frame(X1 = x1, X2 = x2, Y = 0.5+x1+0.6*x2*x2+z+epsilon, Eps=epsilon)
    colnames(df) <-  c("X1","X2","Y", "Eps")
    return(list(df = df, meta = meta, z = z))
  }
  
  y <- Y_true(x1,x2,epsilon)
  df <- data.frame(X1 = x1, X2 = x2, Y = y, Eps = epsilon)
  colnames(df) <-  c("X1","X2","Y", "Eps")
  return(list(df=df,meta=meta))
}

Testowanie poprawności generowania zmiennych o określonych rozkładach

Sprawdźmy działanie naszych funkcji na poszczególnych przypadkach:

Uwaga! Funkcja ggpairs wykonuje tzw. kernel density estimation podczas rysowania funkcji gęstości na podstawie realizacji prób pochodzących z przyjętych teoretycznych rozkładów. Biorąc pod uwagę również to, że nie generujemy bardzo dużej próbki danych dla naszych zmiennych, otrzymujemy przybliżone wykresy gęstości dla rozkładów prawdopodobieństwa ze wstępu na przekątnych widzianych poniżej obrazków.

Wykresy funkcji gęstości zmiennych \(X_1\) i \(X_2\) przedstawiają rozkład normalny ze średnią w punkcie 0. Jednocześnie losowany szum nie zależy od żadnej zmiennej i jego rozkład jest zgodny z przyjętymi założeniami. Generowanie danych jest poprawne dla rozpatrywanego przypadku.

Gęstość prawdopodobieństwa równomiernie rozkłada się dla rozpatrywanych zmiennych objaśniających na przedziale \([-3;3]\). \(\epsilon\) również nie zależy od żadnej zmiennej. Tym samym generowanie danych jest poprawne dla rozpatrywanego przypadku.

Gęstość prawdopodobieństwa równomiernie rozkłada się dla rozpatrywanych zmiennych objaśniających na przedziale \([-5;5]\). \(\epsilon\) również nie zależy od żadnej zmiennej. Tym samym generowanie danych jest poprawne dla rozpatrywanego przypadku.

Wykresy gęstości jądrowej zgadzają się z założonymi rozkładami.

Podobnie jak w poprzednim podpunkcie: generowane dane są zgodne z naszymi założeniami dla tego przypadku. Wartym uwagi jest współczynnik korelacji liniowej Pearsona dla zmiennej \(X_2\) oraz \(Y\), którego wartość jest istotna statystycznie i bliska 1. W rzeczywistości wykres kropkowy przedstawia inną sytuację - zależność jest nieliniowa, a przy większej liczbie obserwacji i przy wzroście wartości zmiennej \(X_2\), \(Y\) wzrasta wykładniczo w stosunku do tej zmiennej - jest to prawdziwa sytuacja.

Z powyższych wykresów widzimy, że rozkład generowanych zmiennych zgadza się z założeniami dla tego przypadku. Zwracamy tu uwagę na dwie mody występujące w rozkładzie zmiennej \(X_1\) oraz na rozkład prawdopodobieństwa dla zmiennej \(X_2\) ze średnią wartością bliską 0.

Widzimy, że wykres gęstości prawdopodobieństwa dla zmiennej \(X_1\) jest zgodny z wykresem rozkładu wykładniczego, a rozkład zmienej \(X_2\) jest rozkładem normalnym ze średnią bliską 0 oraz nieco większym rozrzutem, czyli wariancją danych niż w rozkładzie normalnym standaryzowanym, co właśnie chcieliśmy osiągnąć. Biorąc powyższe pod uwagę oraz to, że korzystaliśmy z funkcji programu R do losowania wartości obserwacji z interesujących nas rozkładów, stwierdzamy, że dane te zostały wygenerowane zgodnie z oczekiwaniami dla tego przypadku.

Widzimy, jak w tym scenariuszu układają się ciężkie ogony przy rozkładzie t-Studenta z trzema stopniami swobody zmiennej \(X_1\), i obserwujemy już wielokrotnie spotykany przedtem przypadek rozkładu normalnego standaryzowanego zmiennej \(X_2\) . Otrzymujemy zatem te rozkłady danych, których się spodziewamy.

Pomijamy opis tego przypadku, ponieważ jest niemalże identyczny do opisu przypadku nr. 6 z jedyną różnicą w tym, że zmienna \(X_2\) pochodzi z rozkładu normalnego, który nie jest standaryzowany. Dane zostały wygenerowany w sposób zgodny z założeniem.

Widzimy charakterystyczne dla rozkładu t-Studenta z trzema stopniami swobody (zmienna \(X_1\)) i dla rozkładu wykładniczego (zmienna \(X_2\) ) wykresy gęstości prawdopodobieństwa. Za każdym razem losowania kolejnych ramek danych otrzymujemy zbliżone wyniki, co sprawia, że dane dla tego przypadku są poprawnie generowane.

Obserwujemy silną dodatnią korelację dla zmiennych \(X_1\) i \(X_2\) pochodzących z dwuwymiarowego rozkładu normalnego. Dane w tym przypadku są generowane zgodnie z oczekiwaniami.

Wykres gęstości prawdopodobieństwa dla zmiennej \(X_2\) wskazuje, że pochodzi ona z rozkładu normalnego standaryzowanego, co się zgadza z naszymi założeniami. Widzimy również statystycznie istotną wartość korelacji pomiędzy zmiennymi \(X_1\) i \(X_2\). Korelacja co prawda jest, ale wiemy, że jest ona nieliniowa (mamy do czynienia z monotoniczną zależnością nieliniową). Nasze założenia są spełnione.

Widzimy niemalże perfekcyjną korelację liniową pomiędzy zmiennymi objaśniającymi, ale wiemy, że kryje się za tym zmienna \(Z\) pochodząca z rozkładu normalnego, która to uzależnia od siebie zmienne \(X_1\) oraz \(X_2\). Szczegóły tego przypadku rozpatrujemy w części 3. tego projektu.

Krótka charakterystyka generowanych danych

Poniższa tabela zawiera empiryczne i teoretyczne informacje dotyczące wygenerowanych ramek danych. Pozwala to na porównanie poprawności otrzymywanych rozkładów.

Empiryczne parametry zmiennych X1 i X2 dla wszystkich rozkładów i powtórzeń
Rozklad mean_X1 expected_mean_x1 var_X1_mean expected_var_x1 mean_X2 expected_mean_x2 var_X2_mean expected_var_x2
double_snorm 0.0116 0.0000 0.9948 1.0000 0.0063 0.0000 0.9918 1.000
double_unif_narrow 0.0064 0.0000 3.0261 3.0000 -0.0045 0.0000 2.9869 3.000
double_unif_wide 0.0021 0.0000 8.2778 8.3330 0.0172 0.0000 8.4556 8.333
norm_exp -0.0253 0.0000 8.7116 8.4100 0.8311 0.8333 0.6993 0.694
snorm_exp 0.0009 0.0000 1.0516 1.0000 0.8262 0.8333 0.6921 0.694
bimodal_snorm 0.4722 0.4200 9.5479 9.4636 -0.0171 0.0000 1.0165 1.000
exp_norm 0.8243 0.8333 0.6496 0.6940 -0.0072 0.0000 8.4423 8.410
t_snorm 0.0332 0.0000 2.8118 3.0000 0.0272 0.0000 1.0145 1.000
t_norm -0.0044 0.0000 2.7619 3.0000 0.0113 0.0000 8.5536 8.410
t_exp -0.0480 0.0000 3.2950 3.0000 0.8366 0.8333 0.7186 0.694
linear_corr_snorm 0.0081 0.0000 1.0332 1.0000 0.0055 0.0000 1.0400 1.000
nonlinear_corr 3.2725 3.0802 116.2893 80.5394 0.0208 0.0000 0.9905 1.000
hidden_impact 11.9319 0.0000 7.3110 7.3000 14.9647 0.0000 1.9656 1.970

Co oznaczają poszczególne kolumny?

Każdy wiersz powyższej tabeli odpowiada jednemu z przypadków symulacyjnych, a każda kolumna przedstawia empiryczne i teoretyczne parametry zmiennych \(X_1\) ​ i \(X_2\) ​ obliczone na podstawie 30 powtórzeń generowania ramek danych z tymi zmiennymi oraz na podstawie teoretycznych obliczeń odpowiednio.

  • Rozklad – nazwa przypadku symulacyjnego, np. t_norm, pois_snorm, bimodal_snorm itd.

  • mean_X1 – średnia wartość empiryczna zmiennej \(X_1\) , obliczona jako średnia ze średnich z 30 wygenerowanych prób

  • expected_mean_x1 – teoretyczna wartość oczekiwana zmiennej \(X_1\)

  • var_X1_mean – średnia z oszacowanych wariancji \(X_1\) z 30 ramek danych

  • expected_var_x1 – teoretyczna wariancja zmiennej \(X_1\)

  • mean_X2 – średnia empiryczna zmiennej \(X_2\) (średnia ze średnich z 30 wygenerowanych ramek danych)

  • expected_mean_x2 – teoretyczna wartość oczekiwana zmiennej \(X_2\)

  • var_X2_mean – średnia z oszacowanych wariancji \(X_2\) z 30 ramek danych

  • expected_var_x2 – teoretyczna wariancja zmiennej \(X_2\)

Przypominamy raz jeszcze, że ta sekcja dotyczyła wyłącznie testowania funkcji do generowania danych. Uzyskane w tej sekcji wyniki dają nam praktyczny wgląd na sposób generowania danych przez dedykowaną funkcję. Widzimy, że rozkłady zgadzają się z założeniami, a tworzone ramki danych są miarodajne. Możemy zatem przejść do dalszych badań.

Modele regresji dla różnych rozkładów danych i ich wpływ na MSE

Dlaczego stosujemy modele liniowe lm()?

Model lm() w R implementuje klasyczną metodę najmniejszych kwadratów (MNK), która znajduje takie współczynniki \(\beta\), aby minimalizować sumę kwadratów błędów między wartościami obserwowanymi (Y) a wartościami przewidywanymi przez model \((\hat{Y})\).
Oznacza to, że lm() szuka prostych (hiperpłaszczyzn) najlepiej dopasowanych do danych w sensie średniokwadratowym (czyli właśnie tego, co mierzy MSE).

Funkcja prawdziwa postaci \(Y=0.5+X_1+0.6*X_2^2 +\epsilon\) jest modelem liniowym względem parametrów, ponieważ zależność od \(\beta_0\), \(\beta_1\), \(\beta_2\) jest liniowa (nawet jeśli występuje kwadrat zmiennej \(X_2^2\)). Dlatego możemy używać lm() nawet mimo obecności nieliniowości w zmiennych.

Dlaczego trzy modele i jakie mają znaczenie?

Wybraliśmy trzy różne modele aby ocenić, jak różny stopień dopasowania modelu do prawdziwej zależności wpływa na wynik metryki MSE.
Każdy z modeli reprezentuje inny poziom trafności założeń:

  • Model Ideal: lm(Y ~ X1 + I(X2^2)) odzwierciedla faktyczną postać funkcji generującej dane (najlepsze możliwe dopasowanie)

  • Model Simple: lm(Y ~ X1 + X2) jest celowo uproszczony, aby pokazać, jak wzrasta błąd, gdy pominiemy ważny nieliniowy element - służy jako punkt odniesienia.

  • Model Mixed: lm(Y ~ X1 + X2 + I(X1*X2)) model bardziej złożony niż klasyczny model liniowy

I() mówi funkcji lm() lub innej modelującej funkcji, że to, co jest w środku, ma być zinterpretowane dosłownie jako wyrażenie matematyczne, a nie jako specjalna składnia formuły R.

Tworzenie modeli i wyznaczanie MSE

Pokaż/Ukryj kod 
library(dplyr)
library(ggplot2)
library(MASS)
library(broom)   
library(gridExtra)

set.seed(123)

cases <- c(
  "double_snorm",
  "double_unif_narrow",
  "double_unif_wide",
  "norm_exp",
  "snorm_exp",
  "bimodal_snorm",
  "exp_norm",
  "t_snorm",
  "t_norm",
  "t_exp",
  "linear_corr_snormrho0.2",
  "linear_corr_snorm_rho0.5",
  "linear_corr_snorm_rho0.9",
  "nonlinear_corr_a0.5b0.2",
  "nonlinear_corr_a1b0.5",
  "nonlinear_corr_a1.5b1",
  "hidden_impact"
)

n_iter <- 200
n_cases <- length(cases)
n_obs <- 120
Pokaż/Ukryj kod 
# Zapis przykładowych ramek danych (po jednej dla przykładu)
if (!dir.exists("data_samples")) dir.create("data_samples")

for (case_name in cases) {
  if (grepl("linear_corr_snorm", case_name)) {
    rho_val <- as.numeric(sub(".*rho", "", case_name))
    df_list <- generate_data(case = "linear_corr_snorm", rho = rho_val)
  } else if (grepl("nonlinear_corr", case_name)) {
    ab_vals <- regmatches(case_name, gregexpr("[0-9\\.]+", case_name))[[1]]
    a_val <- as.numeric(ab_vals[1])
    b_val <- as.numeric(ab_vals[2])
    df_list <- generate_data(case = "nonlinear_corr", a = a_val, b = b_val)
  } else {
    df_list <- generate_data(case = case_name)
  }
  df <- df_list$df
  write.csv(df, paste0("data_samples/", case_name, ".csv"), row.names = FALSE)
}

Podział na zbiór treningowy i testowy

Każdy zbiór danych jest dzielony na:

  • 80% train - do dopasowania modelu (lm())

  • 20% test - do oceny jakości dopasowania (liczymy MSE)

Jeśli policzylibyśmy błąd MSE na tych samych danych, na których model był trenowany, mielibyśmy zaniżony błąd (overfitting).

Podział pozwala realnie ocenić zdolność modelu do generalizacji, czyli jak radzi sobie z nowymi danymi.

W taki sposób podzielone zostały wszystkie ramki danych (200) dla każdego przypadku.

Rola wariancji zmiennych

MSE wyraża rozproszenie danych od linii regresji.

  • Mniejsze MSE oznacza, że model jest lepiej dopasowany, błędy predykcji są mniejsze

  • Większe MSE oznacza, że model gorzej odwzorowuje zależności w danych.

Należy być jednak ostrożnym, ponieważ MSE jest miarą bezwzględną - w zależności od sytuacji znaczenie wielkości może się zmieniać. Przy analizie konkretnych przypadków będziemy zwracać uwagę również na wariancję zmiennych, która wpływa na wartość błędu średniokwadratowego.

MSE nie jest bezpośrednio porównywalne między przypadkami, jeśli skale \(X_{1}\) i \(X_2\) różnią się.

Przykładowo dane o dużym rozrzucie np. \(X\sim N(0,10)\) mogą mieć naturalnie większe błędy prognozy niż dane z \(X \sim N(0,1)\), nawet jeśli dopasowanie względne jest podobne. Dlatego przy analizie należy wziąć pod uwagę wariancję \(X\).

Pokaż/Ukryj kod 
# dopasowanie modeli i liczenie MSE
results_list <- list()
idx <- 1

for (case_id in 1:n_cases) {
  case_name <- cases[case_id]
  
  for (i in 1:n_iter) {
    # generowanie danych
    if (grepl("linear_corr_snorm", case_name)) {
      # wyciągamy rho z nazwy, np. z 'rho0.5'
      rho_val <- as.numeric(sub(".*rho", "", case_name))
      df_list <- generate_data(case = "linear_corr_snorm", rho = rho_val)
      
    } else if (grepl("nonlinear_corr", case_name)) {
      # wyciągamy a i b z nazwy, np. z 'a1.5b1'
      ab_vals <- regmatches(case_name, gregexpr("[0-9\\.]+", case_name))[[1]]
      a_val <- as.numeric(ab_vals[1])
      b_val <- as.numeric(ab_vals[2])
      df_list <- generate_data(case = "nonlinear_corr", a = a_val, b = b_val)
    } else {
      df_list <- generate_data(case = case_name)
    }
    df <- df_list$df
    
    # Split 80/20
    train_idx <- sample(1:nrow(df), size = 0.8 * nrow(df))
    train <- df[train_idx, ]
    test <- df[-train_idx, ]
    
    m_simple <- lm(Y ~ X1 + X2, data = train)
    m_ideal  <- lm(Y ~ X1 + I(X2^2), data = train)
    m_mixed  <- lm(Y ~ X1 + X2 + I(X1 * X2), data = train)
    
    # wyliczanie MSE i zapis współczynników
    models <- list(simple = m_simple, ideal = m_ideal, mixed = m_mixed)
    for (m_name in names(models)) {
      m <- models[[m_name]]      
      preds <- predict(m, newdata=test)
      mse <- mean((test$Y - preds)^2)
      
      coefs <- coef(m)
      row <- tibble(
        case = case_name,
        iter = i,
        model = m_name,
        mse = mse,
        intercept = ifelse("(Intercept)" %in% names(coefs), coefs["(Intercept)"], NA_real_),
        b_X1 = ifelse("X1" %in% names(coefs), coefs["X1"], NA_real_),
        b_X2 = ifelse("X2" %in% names(coefs), coefs["X2"], NA_real_),
        b_X2sq = ifelse("I(X2^2)" %in% names(coefs), coefs["I(X2^2)"], NA_real_),
        b_X1X2 = ifelse("I(X1 * X2)" %in% names(coefs), coefs["I(X1 * X2)"], NA_real_)
      )
      
      results_list[[idx]] <- row
      idx <- idx + 1
    }
  }
}

results_all <- bind_rows(results_list)

# Zapis wyników
write.csv(results_all, "wyniki_parametry_MSE.csv", row.names = FALSE)

Statystyki MSE

MSE jest zmienną losową, ponieważ zależy od losowo wygenerowanych danych. Powtarzamy eksperyment 200 razy dla każdego przypadku, co pozwala nam policzyć

  • średnie - przeciętny błąd

  • mediany - błąd “typowy”, bardziej odporny na wartości odstające

  • kwartyle Q25 i Q75, które pokazują rozrzut błędów

  • odchylenia standardowe - mówi o zmienności wyników między iteracjami

i ocenić stabilność modelu - czy między iteracjami MSE się waha, czy jest stabilne. Sprawdzamy w ten sposób w jakich rozkładach błędy są największe i jak duże znaczenie w dopasowaniu się do danych miał dobór wzoru modelu.

Pokaż/Ukryj kod 
# Statystyki MSE
results_summary <- results_all %>%
  group_by(case, model) %>%
  summarise(
    min_MSE=min(mse),
    q25 = quantile(mse, 0.25),
    mean_MSE = mean(mse),
    median_MSE = median(mse),
    q75 = quantile(mse, 0.75),
    max_MSE=max(mse),
    sd_MSE = sd(mse),
    .groups = "drop"
  )

write.csv(results_summary, "summary_MSE.csv", row.names = FALSE)

results_summary %>%
  mutate(across(where(is.numeric), ~round(., 4))) %>%
  kable(caption = "Statystyki MSE dla poszczególnych przypadków i modeli") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed")) %>%
  column_spec(1, bold = TRUE) %>%
  collapse_rows(columns = 1, valign = "middle")
Statystyki MSE dla poszczególnych przypadków i modeli
case model min_MSE q25 mean_MSE median_MSE q75 max_MSE sd_MSE
bimodal_snorm ideal 0.0042 0.0082 0.0102 0.0096 0.0120 0.0189 0.0029
mixed 0.1354 0.4141 0.7973 0.6696 0.9658 3.7760 0.5559
simple 0.1473 0.4101 0.7559 0.5939 0.9418 4.1909 0.5317
double_snorm ideal 0.0035 0.0083 0.0104 0.0102 0.0126 0.0191 0.0030
mixed 0.1561 0.4176 0.8112 0.6740 1.0470 2.8149 0.5103
simple 0.1596 0.4133 0.7927 0.6619 1.0557 3.1058 0.5136
double_unif_narrow ideal 0.0034 0.0085 0.0107 0.0105 0.0128 0.0188 0.0029
mixed 1.2042 2.3274 2.7298 2.6784 3.0705 4.6891 0.6333
simple 1.2037 2.2756 2.6996 2.6607 3.0154 4.5204 0.6129
double_unif_wide ideal 0.0028 0.0083 0.0105 0.0103 0.0125 0.0221 0.0032
mixed 9.2860 18.1913 21.7211 21.5384 25.1638 39.0761 5.4007
simple 9.8518 18.0734 21.2760 21.6044 24.2744 37.9374 4.9723
exp_norm ideal 0.0047 0.0080 0.0102 0.0101 0.0122 0.0196 0.0029
mixed 9.7640 26.9976 55.1848 45.8961 62.8221 315.2500 46.9501
simple 9.1819 25.8804 52.6158 43.2131 62.5845 314.3302 42.3956
hidden_impact ideal 0.0041 0.0085 0.0105 0.0101 0.0121 0.0214 0.0031
mixed 0.0041 0.0085 0.0105 0.0101 0.0121 0.0214 0.0031
simple 0.4090 1.5255 3.1123 2.4248 3.8655 14.7734 2.2301
linear_corr_snorm_rho0.5 ideal 0.0043 0.0080 0.0102 0.0098 0.0122 0.0201 0.0031
mixed 0.1162 0.2725 0.4740 0.3749 0.5575 2.0048 0.3348
simple 0.1394 0.3853 0.7723 0.5959 0.9238 3.7862 0.5961
linear_corr_snorm_rho0.9 ideal 0.0044 0.0083 0.0105 0.0105 0.0128 0.0192 0.0029
mixed 0.0215 0.0587 0.0904 0.0815 0.1144 0.2518 0.0448
simple 0.1730 0.3888 0.7513 0.6216 0.9944 2.6078 0.4840
linear_corr_snormrho0.2 ideal 0.0043 0.0084 0.0104 0.0103 0.0122 0.0187 0.0027
mixed 0.1561 0.4055 0.7606 0.5901 0.9027 3.1477 0.5390
simple 0.1749 0.4532 0.7926 0.6427 0.9488 3.3684 0.5357
nonlinear_corr_a0.5b0.2 ideal 0.0036 0.0080 0.0102 0.0098 0.0121 0.0198 0.0029
mixed 0.0552 0.1596 0.2614 0.2192 0.3074 1.4891 0.1718
simple 0.1673 0.4194 0.7529 0.6097 0.8760 4.5577 0.5439
nonlinear_corr_a1.5b1 ideal 0.0042 0.0081 0.0102 0.0098 0.0121 0.0228 0.0031
mixed 0.0202 0.0415 0.0751 0.0563 0.0749 0.8991 0.0802
simple 0.0200 0.0911 0.2847 0.1427 0.2493 4.0425 0.5114
nonlinear_corr_a1b0.5 ideal 0.0047 0.0086 0.0104 0.0097 0.0121 0.0187 0.0028
mixed 0.0154 0.0553 0.1458 0.0937 0.1452 2.1437 0.2036
simple 0.0601 0.1383 0.2110 0.1764 0.2370 1.9686 0.1618
norm_exp ideal 0.0037 0.0078 0.0102 0.0100 0.0121 0.0272 0.0034
mixed 0.0537 0.1902 0.8671 0.3668 0.6923 19.8119 2.0971
simple 0.0522 0.1951 0.7535 0.3549 0.6682 15.8515 1.6865
snorm_exp ideal 0.0045 0.0083 0.0100 0.0098 0.0114 0.0188 0.0026
mixed 0.0441 0.2018 0.7617 0.3878 0.8211 15.0120 1.3004
simple 0.0508 0.2041 0.6952 0.3686 0.7183 14.9286 1.2711
t_exp ideal 0.0042 0.0083 0.0103 0.0099 0.0120 0.0231 0.0031
mixed 0.0358 0.1933 0.7128 0.3281 0.5861 11.2516 1.3474
simple 0.0353 0.1931 0.6470 0.3241 0.5666 10.9228 1.2003
t_norm ideal 0.0037 0.0086 0.0107 0.0102 0.0126 0.0209 0.0031
mixed 10.6921 32.0753 58.9187 47.1701 67.8167 589.9308 53.0344
simple 10.6071 28.9494 57.0330 46.4800 68.3398 593.1615 52.4424
t_snorm ideal 0.0042 0.0082 0.0102 0.0101 0.0121 0.0185 0.0029
mixed 0.1630 0.4268 0.7944 0.6313 0.9302 4.2049 0.5745
simple 0.1562 0.4159 0.7559 0.6070 0.9392 3.0123 0.4905

Zwróćmy uwagę na średnią wartość MSE. Dla wszystkich przypadków model Ideal czyli lm(Y~X1+I(X2^2)) osiągał najniższe średnie wartości błędu MSE i najmniejszą zmienność błędu, co jest naturalne, ponieważ odpowiada on dokłądnie generującej funkcji. Błąd na poziomie 0.01 jest głównie spowodowany szumem losowym więc jest to najlepsze możliwe dopasowanie. W rzeczywistości jednak, raczej nie znamy funkcji generującej dane, a szukamy modelu, który najlepiej dopasuje się do danych. Tak traktujemy dwa pozostałe modele Mixed lm(Y~X1+X2+I(X1*X2)) i Simple lm(Y~X1+X2) . Dla tych modeli obserwujemy większe zróżnicowanie wartości MSE - między różnymi rozkładami, ale również w obrębie jednego przykładu. Przyjrzyjmy się kilku przypadkom, aby lepiej zrozumieć skąd te różnice się biorą.

Pokaż/Ukryj kod 
# Wykres: średnie MSE (słupki)
ggplot(results_summary, aes(x = case, y = mean_MSE, fill = model)) +
  geom_col(position = "dodge") +
  theme_minimal() +
  labs(title = "Średnie wartości MSE w poszczególnych przypadkach",
       x = "Przypadek",
       y = "Średnie MSE") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Z wykresu widzimy, że dla 4 przypadków rozkładów modele mixed i simple poradziły sobie bardzo źle. Dla porównania zestawimy wyniki dla przypadków:

  • double_unif_narrow i double_unif_wide

  • t_norm i t_snorm

  • exp_norm i norm_exp

W tym celu wczytane zostaną dane z 1 ramki dla każdego przykładu

Pokaż/Ukryj kod 
plot_case_fit <- function(case_name) {
  df <- read.csv(paste0("data_samples/", case_name, ".csv"))
  
  set.seed(123)
  train_idx <- sample(1:nrow(df), size = 0.8 * nrow(df))
  train <- df[train_idx, ]
  test <- df[-train_idx, ]
  
  m_simple <- lm(Y ~ X1 + X2, data = train)
  m_ideal  <- lm(Y ~ X1 + I(X2^2), data = train)
  m_mixed  <- lm(Y ~ X1 + X2 + I(X1 * X2), data = train)
  
  test$pred_simple <- predict(m_simple, newdata = test)
  test$pred_ideal  <- predict(m_ideal, newdata = test)
  test$pred_mixed  <- predict(m_mixed, newdata = test)
  
  mse_simple <- mean((test$Y - test$pred_simple)^2)
  mse_ideal  <- mean((test$Y - test$pred_ideal)^2)
  mse_mixed  <- mean((test$Y - test$pred_mixed)^2)
  var_X1 <- var(test$X1)
  var_X2 <- var(test$X2)
  var_Y  <- var(test$Y)
  
  summary_table <- tibble(
    case = case_name,
    var_X1 = var_X1,
    var_X2 = var_X2,
    var_Y = var_Y,
    mse_simple = mse_simple,
    mse_ideal = mse_ideal,
    mse_mixed = mse_mixed
  )
  
  # Wykresy
  p1 <- ggplot(test, aes(x = X1, y = Y)) +
    geom_point(alpha = 0.5, color = "grey40") +
    geom_line(aes(y = pred_simple, color = "Simple"), linewidth = 0.8) +
    geom_line(aes(y = pred_ideal, color = "Ideal"), linewidth = 0.8) +
    geom_line(aes(y = pred_mixed, color = "Mixed"), linewidth = 0.8) +
    labs(
      title = paste("Dopasowanie względem X1"),
      x = "X1", y = "Y", color = "Model"
    ) +
    theme_minimal()
  
  p2 <- ggplot(test, aes(x = X2, y = Y)) +
    geom_point(alpha = 0.5, color = "grey40") +
    geom_line(aes(y = pred_simple, color = "Simple"), linewidth = 0.8) +
    geom_line(aes(y = pred_ideal, color = "Ideal"), linewidth = 0.8) +
    geom_line(aes(y = pred_mixed, color = "Mixed"), linewidth = 0.8) +
    labs(
      title = paste("Dopasowanie względem X2"),
      x = "X2", y = "Y", color = "Model"
    ) +
    theme_minimal()
  
  print(case_name)
  
  grid.arrange(p1, p2, ncol = 2)
  
  summary_table %>%
    kbl(caption = paste("Statystyki dla przypadku:", case_name)) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  full_width = FALSE)
}
Pokaż/Ukryj kod 
plot_pred_vs_actual <- function(case_name) {
  df <- read.csv(paste0("data_samples/", case_name, ".csv"))
  
  set.seed(123)
  train_idx <- sample(1:nrow(df), size = 0.8 * nrow(df))
  train <- df[train_idx, ]
  test <- df[-train_idx, ]
  
  m_simple <- lm(Y ~ X1 + X2, data = train)
  m_ideal  <- lm(Y ~ X1 + I(X2^2), data = train)
  m_mixed  <- lm(Y ~ X1 + X2 + I(X1 * X2), data = train)
  
  # Predykcje
  test$pred_simple <- predict(m_simple, newdata = test)
  test$pred_ideal  <- predict(m_ideal, newdata = test)
  test$pred_mixed  <- predict(m_mixed, newdata = test)
  
  # Obliczanie MSE
  mse_simple <- mean((test$Y - test$pred_simple)^2)
  mse_ideal  <- mean((test$Y - test$pred_ideal)^2)
  mse_mixed  <- mean((test$Y - test$pred_mixed)^2)
  
  # Pomocnicza funkcja do tworzenia wykresu dla jednego modelu
  make_plot <- function(pred_col, model_name, mse_value, color) {
    ggplot(test, aes(x = Y, y = !!sym(pred_col))) +
      geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "black") +
      geom_segment(aes(xend = Y, yend = !!sym(pred_col)), alpha = 0.4, color = color) +
      geom_point(color = color, size = 1.5) +
      labs(
        title = paste("Model", model_name),
        subtitle = sprintf("MSE = %.4f", mse_value),
        x = "Prawdziwy Y",
        y = "Przewidziany Ŷ"
      ) +
      theme_minimal() +
      coord_equal()
  }
  
  # Trzy wykresy obok siebie
  p1 <- make_plot("pred_simple", "Simple", mse_simple, "darkred")
  p2 <- make_plot("pred_ideal", "Ideal", mse_ideal, "steelblue")
  p3 <- make_plot("pred_mixed", "Mixed", mse_mixed, "darkgreen")
  
  print(case_name)
  
  grid.arrange(p1, p2, p3, nrow = 1)
}

Double Uniform Narrow vs. Double Uniform Wide

Obie wersje mają ten sam rodzaj rozkładu - jednostajny - ale różny zakres:

  • narrow: \(X_i \sim U(-3,3)\)

  • wide: \(X_i \sim U(-5,5)\)

Pokaż/Ukryj kod 
plot_case_fit("double_unif_narrow")
[1] "double_unif_narrow"

Statystyki dla przypadku: double_unif_narrow
case var_X1 var_X2 var_Y mse_simple mse_ideal mse_mixed
double_unif_narrow 3.172489 3.013163 5.848907 3.165899 0.0081483 2.913976
Pokaż/Ukryj kod 
plot_pred_vs_actual("double_unif_narrow")
[1] "double_unif_narrow"

Na podstawie powyższego wykresu scatterplot możemy powiedzieć jak dobrze model dopasowuje się do danych. Im bliżej linii \(\hat{Y} = Y\) znajdują się punkty przedstawiające wartości przewidywane, tym lepszy jest model. Widać wyraźnie, że model Ideal “rzuca” przewidywania w większości pokrywające się z prawdziwymi wartościami.

Pokaż/Ukryj kod 
diagnose_model <- function(case_name) {
  df <- read.csv(paste0("data_samples/", case_name, ".csv"))
  
  model <- lm(Y ~ X1 + I(X2^2), data = df)
  
  par(mfrow = c(2, 2))
  plot(model, main = paste("Diagnostyka modelu:", case_name))
  par(mfrow = c(1, 1)) 
}
Pokaż/Ukryj kod 
diagnose_model("double_unif_narrow")

Dla idealnie dopasowanego modelu wykres Residuals vs Fitted (reszty względem wartości dopasowanych) nie powinien wykazywać żadnego widocznego wzoru. Punkty powinny być rozmieszczone losowo wokół osi poziomej. Jeśli natomiast zauważylibyśmy np. kształt przypominający parabolę, mogłoby to oznaczać, że zależność między zmiennymi nie jest liniowa.
W naszym przypadku punkty są rozproszone nieregularnie, co sugeruje, że przyjęty model liniowy jest odpowiedni.

Wykres Normal Q-Q służy do oceny, czy reszty (błędy) mają rozkład normalny. Jeśli pochodzą z rozkładu normalnego, punkty powinny układać się wzdłuż prostej przerywanej.
W rozważanym przypadku obserwujemy, że punkty przebiegają wzdłuż tej linii, co potwierdza, że błędy mają rozkład zbliżony do normalnego.

Wykres Scale-Location (znany też jako Spread-Location) pokazuje, czy wariancja błędów jest stała (homoscedastyczność). Jeśli czerwona linia trendu jest prawie pozioma, a punkty są równomiernie rozproszone, możemy uznać, że założenie jednorodności wariancji jest spełnione.

Wykres Residuals vs Leverage pomaga zidentyfikować obserwacje wpływowe, czyli takie, które znacząco oddziałują na estymowane parametry modelu. Punkty znajdujące się poza liniami Cook’s distance (zaznaczonymi przerywanymi liniami) są potencjalnie wpływowe - ich usunięcie mogłoby istotnie zmienić dopasowanie modelu.
Nasz model nie posiada potencjalnie wpływowych obserwacji.

Pokaż/Ukryj kod 
plot_case_fit("double_unif_wide")
[1] "double_unif_wide"

Statystyki dla przypadku: double_unif_wide
case var_X1 var_X2 var_Y mse_simple mse_ideal mse_mixed
double_unif_wide 8.672501 9.898818 37.90714 24.04106 0.0102462 28.77621
Pokaż/Ukryj kod 
plot_pred_vs_actual("double_unif_wide")
[1] "double_unif_wide"

Pokaż/Ukryj kod 
diagnose_model("double_unif_wide")

Z powyższych wykresów możemy wyciągnąć podobne wnioski co w przypadku double_unif_narrow.

“double_unif” model narrow wide
Simple 3.121 23.947
Ideal 0.011 0.009
Mixed 2.858 28.41

Model idealny w obu przypadkach ma bardzo niskie MSE, a to oznacza, że poprawnie odwzorowuje funkcję generującą dane niezależnie od zakresu zmiennych.

Natomiast modele Simple i Mixed mają MSE znacznie większe dla „wide”, ponieważ wartości \(X_2\) i \(X_1\) są rozrzucone szerzej, co skutkuje większym zróżnicowaniem wartości \(Y\), a przez to większymi błędami predykcji - bo błędy kwadratowe \((Y - \hat{Y})^2\) rosną szybciej, im większy rozrzut (wariancja \(X_1\), \(X_2\)). Wartości \(Y\) też przyjmują większy zakres, bo są funkcją \(X_1\) i \(X_2\).

Wariancje zmiennych \(X_1\) i \(X_2\) w przypadkach narrow i wide różnią się kilkukrotnie. To oznacza, że skala błędów w modelach prostych (Simple, Mixed) rośnie proporcjonalnie do rozrzutu zmiennych, ponieważ MSE jest miarą bezwzględną, a nie znormalizowaną.

Sama skala i wariancja zmiennych wpływa na wartość MSE, nawet jeśli dopasowanie względne (np. \(R^2\)) jest podobne. Dlatego MSE nie jest bezwzględnie porównywalne między przypadkami o różnych skalach danych.

(Współczynnik determinacji \(R^2\) pokazuje, jak dobrze model wyjaśnia zmienność danych niezależnie od ich skali. MSE natomiast rośnie wraz ze skalą i wariancją zmiennych, dlatego przy porównywaniu modeli dla różnych rozkładów bardziej miarodajne jest porównanie \(R^2\) niż samych wartości MSE.

Przykład z rozkładem t-Studenta (cięższe ogony)

Porównujemy t_norm i t_snorm, dla których w obu przypadkach \(X_1 \sim t(3)\) czyli są to rozkłady t-Studenta z 3 stopniami swobody. Natomiast wariancja \(X_2\) się różni. t_snorm \(X_2 \sim N(0,1)\), a t_snorm \(X_2 \sim N(0,2.9)\)

Dla t_snorm obie zmienne mają umiarkowany rozrzut. Zmienna \(X_2\) w t_norm ma znacznie większe odchylenie standardowe, co oznacza, że dane są bardziej rozproszone w przestrzeni.

Pokaż/Ukryj kod 
plot_case_fit("t_norm")
[1] "t_norm"

Statystyki dla przypadku: t_norm
case var_X1 var_X2 var_Y mse_simple mse_ideal mse_mixed
t_norm 2.310711 6.934366 20.48763 19.13618 0.0081343 23.31183
Pokaż/Ukryj kod 
plot_pred_vs_actual("t_norm")
[1] "t_norm"

Pokaż/Ukryj kod 
diagnose_model("t_norm")

Pokaż/Ukryj kod 
plot_case_fit("t_snorm")
[1] "t_snorm"

Statystyki dla przypadku: t_snorm
case var_X1 var_X2 var_Y mse_simple mse_ideal mse_mixed
t_snorm 1.183799 0.7243344 1.341954 0.2843283 0.0112532 0.3059577

Dla t_snorm wartości \(X_2\) są małe (blisko zera), więc składnik \(X_2^2\) nie ma dużego wpływu na \(Y\) -> Model łatwo się dopasowuje, błędy są małe.

Dla t_norm, wartości \(X_2\) ​są dużo większe - kwadrat rośnie szybko więc różnice w \(Y\) są dużo większe -> nawet mały błąd prognozy \(\hat{Y}\) przekłada się na duży błąd kwadratowy, dlatego MSE rośnie wykładniczo wraz z rozrzutem zmiennych.

model t_snorm t_norm
Simple 0.446 141.746
Ideal 0.01 0.01
Mixed 0.46 147.152
Pokaż/Ukryj kod 
plot_pred_vs_actual("t_snorm")
[1] "t_snorm"

Pokaż/Ukryj kod 
diagnose_model("t_snorm")

W t_norm zobaczymy większy rozrzut reszt, co wskazuje na heteroskedastyczność (rosnącą wariancję błędów przy większych wartościach X2).

Wariant t_norm ma znacznie większy rozrzut - większą wariancję \(X_2\) i cięższe ogony \(X_1\). Oznacza to, że w danych częściej pojawiają się wartości skrajne (outliery), przez co błędy kwadratowe przy modelach niedopasowanych (Simple, Mixed) są ekstremalnie duże.

Dla modelu Ideal różnice są niewielkie, bo model dokładnie odzwierciedla funkcję generującą dane. Dla modeli Simple i Mixed MSE dla t_norm jest kilkaset razy większe, mimo że dopasowanie względne (np. \(R^2\) ) może być podobne.

To dowód, że MSE jest wrażliwe na skalę danych i nie nadaje się do bezpośredniego porównywania między przypadkami o różnych wariancjach.

Wniosek: Rozkłady z większą wariancją lub cięższymi ogonami (np. t-Studenta, szeroki uniform) powodują znaczne zwiększenie MSE przy modelach uproszczonych, mimo identycznego modelu funkcji.

Norm_exp i exp_norm

Przy norm_exp \(X_1 \sim N(0, 2.9)\) i \(X_2 \sim Exp(1.2)\). Więc jedna zmienna jest symetryczna, a druga skośna w prawo.

Pokaż/Ukryj kod 
plot_case_fit("norm_exp")
[1] "norm_exp"

Statystyki dla przypadku: norm_exp
case var_X1 var_X2 var_Y mse_simple mse_ideal mse_mixed
norm_exp 7.3134 1.201378 7.41053 0.2582772 0.0078207 0.2438117
Pokaż/Ukryj kod 
plot_pred_vs_actual("norm_exp")
[1] "norm_exp"

Oczywiście obserwujemy bardzo dobre dopasowanie modelu Ideal, ale co ciekawe modele Mixed i Simple mają małe MSE. Wynika to z zakresu i rozkładu \(X_2\):

\(X_2\) z rozkładu wykładniczego o \(\lambda =1.2\) przyjmuje małe wartości (średnio ok. 0.83), a większe wartości są bardzo rzadkie. W tym zakresie \(X_2^2 \approx X_2\) (bo przy małych wartościach, kwadrat nie zmienia dużo.

Czyli model \(Y \sim X_1 +X_2\) nie „widzi” mocno nieliniowości - i mimo że teoretycznie jest źle określony, praktycznie dopasowuje się bardzo dobrze, bo dane nie pokazują dużych odchyleń od liniowości.

Model \(Y \sim X_1 +X_2+X1\cdot X2\) osiągnął wartość MSE bardzo podobną do Simple. Składnik \(X1\cdot X2\) jest niepotrzebny, jednak nie ma on dużego wpływu na model, jeśli jego współczynnik jest mały.

Pokaż/Ukryj kod 
diagnose_model("norm_exp")

  • Residuals vs Fitted:
    Widoczny lekki wzór - reszty nie są całkowicie losowe, co sugeruje,
    że model liniowy może nie do końca uchwycić nieliniową zależność generowaną przez \(X_2^2\).
  • Scale-Location:
    Czerwona linia jest pozioma i punkty równomiernie rozrzucone - wariancje są jednorodne.
  • Residuals vs Leverage:
    Pojedyncze punkty blisko granicy - czyli dane z długiego ogona rozkładu wykładniczego
    mają większy wpływ na dopasowanie.

W przypadku exp_norm \(X_1 \sim Exp(1.2)\) i \(X_2 \sim N(0, 2.9)\)

Pokaż/Ukryj kod 
plot_case_fit("exp_norm")
[1] "exp_norm"

Statystyki dla przypadku: exp_norm
case var_X1 var_X2 var_Y mse_simple mse_ideal mse_mixed
exp_norm 1.114356 9.109857 107.2702 108.5138 0.0079043 124.1917
Pokaż/Ukryj kod 
plot_pred_vs_actual("exp_norm")
[1] "exp_norm"

Gdy zamienimy rozkłady, z których pochodzą zmienne objaśniające, drastycznie zwiększa się wartość błędu średniokwadratowego dla modeli Simple i Mixed.

W tym wypadku składnik \(0.6\cdot X_2^2\) ma szeroki rozrzut. Składnik \(0.6\cdot X_2^2\) dominuje w Y, czego prosty model liniowy nie jest w stanie uchwycić.

Pokaż/Ukryj kod 
diagnose_model("exp_norm")

Model nie w pełni dopasował się do liniowej zależności w danych - punkty na 1. wykresie są skupione w okolicach 0.

Dla dużych wartości przewidywanych model ma większy rozrzut błędów -> heteroscedastyczność. Model lepiej radzi sobie dla małych Y, gorzej dla dużych.

Większość punktów ma bardzo małą dźwignię, model nie posiada obserwacji wpływowych.

Oznacza to, że zamiana rozkładów zmiennych objaśniających może diametralnie zmienić charakter zależności i dopasowanie modeli, nawet jeśli sama struktura wzoru pozostaje taka sama.

Zmienne ukryte - analia wpływu na model

Intuicja - przykład z Paryża

Pani Profesor Julie Josse na zimowej szkole z Causality i XAI w Srobonne University, przedstawiła problem zmiennych ukrytych w bardzo intuicyjny i zabawny sposób aby zrozumieć dlaczego jest to istotne i czemu same mlowe i statystyczne podejście może nas zawieść. Możemy np. dowieść statystycznie, że spanie w butach powoduje ból głowy następnego dnia. Pani Profesor dowiodła to testem statystycznym z danych, że ta korelacja jest istotna. W danych nie mieliśmy jednak informacji o tym, że dana osoba spożyła wcześniej spore ilości alkoholu, który sprawił, że po powrocie do domu z imprezy nie była w stanie zdjąć butów i poszła spać. Rano obudziła się z kacem - tego nam same dane nie powiedzą i potrzebujemy causal approach. Posłużyłem się chatem gpt do wygenerowania poniższej grafiki w celu ilustracji, abyśmy zrozumieli lepiej graf przyczynowo skutkowo korelacyjny (chat nie mógł zrozumieć prawidłowo kierunku strzałek więc zostawiam grafikę tak jak ja wygenerował :D).

Przechodząc do konkretnego przypadku naszego projektu w tym podrozdziale badamy przypadek, w którym nieobserwowana zmienna \(Z\) oddziałuje zarówno na zmienne \(X1, X2\)jak i na zmienną objaśnianą \(Y\) Taki układ powoduje:

  • zależność pomiędzy \(X1\) i \(X2\)

  • systematyczny składnik wpływający na \(Y\) którego nie potrafimy „wyjaśnić” przez posiadane zmienne z macierzy eksperymentu, tzw. bias przez pominięcie

  • potencjalnie złudną poprawę na danych uczących i gorszą generalizację.

  • Ponieważ \(Z\) nie jest w rozpatrywanej ramce danych, pokażemy, jak jego wpływ widoczny jest pośrednio: w strukturze predyktorów oraz w resztach modeli oraz w zmianie MSE.

Zbadajmy to jaki “cień” zostawiła w naszych zmiennych \(X-owych\) nieznana zmienna \(Z\) za pomocą techniki PCA (analizy głównych składowych) a jak wiemy główną składową każdej z zmiennych \(X-owych\) jest ukryta i bardzo tajemnicza zmienna \(Z\)

Czy widać tak naprawdę wpływ zmiennej \(Z\) na zmienne \(X-owe\)?

Pokaż/Ukryj kod 
library(dplyr)
library(tibble)    
library(purrr)     
library(ggplot2)    
library(knitr)
library(kableExtra)
library(scales)

dfs <- list(
  df1 = df1_hidden_impact,
  df2 = df2_hidden_impact,
  df3 = df3_hidden_impact
)
Pokaż/Ukryj kod 
pc1_score <- function(data){ pr <- prcomp(scale(dplyr::select(data, X1, X2))) 
list(score = as.numeric(pr$x[,1]), var_expl = pr$sdev[1]^2 / sum(pr$sdev^2)) }


proxy_tbl <- imap_dfr(dfs, function(dfi, nm){
pc1 <- pc1_score(dfi)
tibble(
zbior = nm,
cor_X1_X2 = cor(dfi$X1, dfi$X2),
var_expl_PC1 = round(pc1$var_expl, 3),
cor_PC1_X1 = cor(pc1$score, dfi$X1),
cor_PC1_X2 = cor(pc1$score, dfi$X2)
)
})



proxy_tbl %>%
  transmute(
    `Zbiór`                = zbior,
    `Corr(X1, X2)`         = round(cor_X1_X2, 3),
    `Corr(PC1, X1)`        = round(cor_PC1_X1, 3),
    `Corr(PC1, X2)`        = round(cor_PC1_X2, 3),
    `PC1: udział wariancji`= scales::percent(var_expl_PC1, accuracy = 0.1)
  ) %>%
  kable(align = "lcccc",
        booktabs = TRUE,
        caption = "Sygnały wspólnej przyczyny (proxy-Z)") %>%
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover")) %>%
  column_spec(1, bold = TRUE) %>%
  add_header_above(c(" " = 1, "Korelacje" = 3, " " = 1))
Sygnały wspólnej przyczyny (proxy-Z)
Korelacje
Zbiór Corr(X1, X2) Corr(PC1, X1) Corr(PC1, X2) PC1: udział wariancji
df1 1 -1 -1 100.0%
df2 1 1 1 100.0%
df3 1 -1 -1 100.0%

Oczywiście że tak ponieważ w wariancie hidden_impact przyjęliśmy deterministyczne zależności. Skutkiem jest rząd macierzy cech równy 1, tj. \(X1\) i \(X2\) leżą na jednej prostej, co daje bardzo silne zależności, PC1 odtwarza kierunek \(Z\) z dokładnością do znaku. Jest to zamierzony, skrajny przykład confoundingu: cała zmienność predyktorów pochodzi ze wspólnej, nieobserwowanej przyczyny.

Pokaż/Ukryj kod 
pc1_all <- purrr::imap_dfr(dfs, function(dfi, nm){
  pc <- pc1_score(dfi)
  tibble(
    zbior    = nm,
    PC1      = pc$score,
    var_expl = pc$var_expl
  )
})


facet_labels <- pc1_all %>%
  distinct(zbior, var_expl) %>%
  mutate(label = paste0(zbior, " - PC1: ", percent(var_expl, accuracy = 0.1))) %>%
  { setNames(.$label, .$zbior) }


pc1_hist <- ggplot(pc1_all, aes(x = PC1)) +
  geom_histogram(aes(y = ..density..), bins = 30, alpha = 0.75) +
  geom_density(adjust = 1.2, linewidth = 0.8) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  facet_wrap(~ zbior, scales = "free_x", labeller = as_labeller(facet_labels)) +
  labs(
    title    = "Rozkład PC1 (proxy-Z) w zestawach hidden_impact",
    x        = "PC1 (standaryzowany score z PCA na {X1, X2})",
    y        = "Gęstość"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    strip.text = element_text(face = "bold"),
    plot.title = element_text(face = "bold")
  )

pc1_hist

Wniosek z wykresu jest widoczny gołym okiem ponieważ 100% udziału wskazuje że cała zmienność zmiennych \(X1\) i \(X2\) leży w jednym kierunku. Co potwierdza wniosek z poprzedniej powyższej tabeli.

Dlaczego to ważne w kontekście causal approach?

Gdy PC1 wyjaśnia bardzo dużo wariancji, mamy silne przesłanki, że w danych działa nieobserwowany czynnik.

W takiej sytuacji modele regresji bez tego czynnika będą zostawiać strukturę w resztach (co pokazujemy w kolejnych wykresach) i będą miały podwyższony błąd generalizacji.

To uzasadnia, dlaczego w praktyce warto szukać proxy (zmiennych zastępczych) lub instrumentów żeby zredukować wpływ pominiętej zmiennej.

Walidcja Krzyżowa modeli z ukrytą zmienną wpływową

Ustalmy dwie formy modelu którymi będziemy się posługiwać w dalszych rozważaniach

\[ Y \sim X1 + X2 \]

\[ Y \sim X1 + I(X_2^2) \]

Pokaż/Ukryj kod 
form_simple <- Y ~ X1 + X2
form_ideal <- Y ~ X1 + I(X2^2)
Pokaż/Ukryj kod 
cv_mse <- function(form, data, K = 5, seed = 123){
  set.seed(seed)
  n <- nrow(data)
  folds <- sample(rep(1:K, length.out = n))  # losowe przypisanie do K foldów
  mse_folds <- sapply(1:K, function(k){
    tr <- data[folds != k, , drop = FALSE]
    te <- data[folds == k, , drop = FALSE]
    fit <- lm(form, data = tr)
    mean((te$Y - predict(fit, newdata = te))^2)
  })
  list(mean = mean(mse_folds),
       sd   = sd(mse_folds),
       se   = sd(mse_folds) / sqrt(K),
       per_fold = mse_folds)
}
Pokaż/Ukryj kod 
K <- 5
seed_cv <- 123
cv_tbl <- imap_dfr(dfs, function(dfi, nm){
  cv_s <- cv_mse(form_simple, dfi, K = K, seed = seed_cv)
  cv_i <- cv_mse(form_ideal,  dfi, K = K, seed = seed_cv)
  tibble(
    zbior  = nm,
    model  = c("simple", "ideal"),
    CV_MSE = c(cv_s$mean, cv_i$mean),
    SE     = c(cv_s$se,   cv_i$se)
  )
})

cv_folds_long <- imap_dfr(dfs, function(dfi, nm){
  list(
    tibble(zbior = nm, model = "simple", fold = 1:K,
           mse = cv_mse(form_simple, dfi, K = K, seed = seed_cv)$per_fold),
    tibble(zbior = nm, model = "ideal",  fold = 1:K,
           mse = cv_mse(form_ideal,  dfi, K = K, seed = seed_cv)$per_fold)
  ) |> bind_rows()
})
Pokaż/Ukryj kod 
library(dplyr)
library(tidyr)
library(knitr)
library(kableExtra)
library(scales)

table_pretty <- cv_tbl %>%
  dplyr::select(zbior, model, CV_MSE, SE) %>%
  dplyr::mutate(CI_low = CV_MSE - 1.96 * SE,
                CI_hi  = CV_MSE + 1.96 * SE) %>%
  tidyr::pivot_wider(names_from = model, values_from = c(CV_MSE, SE, CI_low, CI_hi)) %>%
  dplyr::mutate(
    `Δ (simple - ideal)`           = CV_MSE_simple - CV_MSE_ideal,
    `Rel. poprawa ideal vs simple` = ifelse(
      is.finite(CV_MSE_simple) & CV_MSE_simple != 0,
      (CV_MSE_simple - CV_MSE_ideal) / CV_MSE_simple, NA_real_
    )
  )

table_pretty %>%
  dplyr::mutate(
    dplyr::across(
      .cols = c(CV_MSE_simple, SE_simple, CI_low_simple, CI_hi_simple,
                CV_MSE_ideal,  SE_ideal,  CI_low_ideal,  CI_hi_ideal,
                `Δ (simple - ideal)`),
      ~ round(.x, 4)
    ),
    `Rel. poprawa ideal vs simple` = percent(`Rel. poprawa ideal vs simple`, accuracy = 0.1)
  ) %>%
  dplyr::rename(
    Zbiór       = zbior,
    `CV-MSE`    = CV_MSE_simple, `SE`    = SE_simple,
    `CI (low)`  = CI_low_simple, `CI (hi)` = CI_hi_simple,
    `CV-MSE `   = CV_MSE_ideal,  `SE `   = SE_ideal,
    `CI (low) ` = CI_low_ideal,  `CI (hi) ` = CI_hi_ideal
  ) %>%
  knitr::kable(align = "lrrrrrrrrrr", booktabs = TRUE,
               caption = "K-fold CV-MSE (K=5): porównanie modeli bez Z (średnia, SE, 95% CI, różnice)") %>%
  kableExtra::kable_styling(full_width = FALSE, bootstrap_options = c("striped","hover")) %>%
  kableExtra::add_header_above(c(" " = 1, "Model: simple" = 4, "Model: ideal" = 4, "Różnice" = 2)) %>%
  kableExtra::column_spec(1, bold = TRUE)
K-fold CV-MSE (K=5): porównanie modeli bez Z (średnia, SE, 95% CI, różnice)
Model: simple
Model: ideal
Różnice
Zbiór CV-MSE CV-MSE SE SE CI (low) CI (low) CI (hi) CI (hi) Δ (simple - ideal) Rel. poprawa ideal vs simple
df1 4.5094 0.0126 1.4715 0.0020 1.6253 0.0086 7.3935 0.0165 4.4968 99.7%
df2 2.1652 0.0095 0.5558 0.0009 1.0759 0.0077 3.2545 0.0113 2.1557 99.6%
df3 2.5123 0.0105 1.0422 0.0016 0.4695 0.0074 4.5551 0.0135 2.5018 99.6%

Gdzie:

CV_MSE_model - średni błąd walidacji krzyżowej (K-fold) danego modelu.

SE_model – błąd standardowy średniej CV-MSE (zmienność między foldami).

CI_low_model, CI_hi_model – dolna i górna granica 95% CI dla CV-MSE (≈ CV_MSE ± 1.96·SE).

Co z tego wynika dla causal approach?

  • Specyfikacja modelu ma pierwszorzędne znaczenie. Sama obecność nieobserwowanej \(Z\) nie musi niszczyć jakości predykcji, jeśli informacja o \(Z\) jest „pośrednio zakodowana” w \(X\) i uchwycona właściwą postacią

  • Brak właściwej postaci ma wpływ duży koszt generalizacji. simple traci nie dlatego, że nie widzi \(Z\), tylko dlatego, że źle modeluje zależność \(Y\) do \(X\).

Podsumowanie sekcji zmiennych ukrytych

Konkluzja. W naszym skrajnym scenariuszu hidden_impact dominująca struktura w \(X\) (PC1 ≈ cień \(Z\)) oraz CV-MSE pokazują, że kluczowa jest prawidłowa postać funkcjonalna względem obserwowalnych cech; dopiero gdy po jej uwzględnieniu reszty nadal są skorelowane z proxy-Z i CV-MSE pozostaje wysokie, sensowne jest inwestowanie w dodatkowe zmienne/proxy/instrumenty.

Pokaż/Ukryj kod 
plot_case_fit2 <- function(case_name) {
  df <- read.csv(paste0("data_samples/", case_name, ".csv"))
  
  set.seed(123)
  train_idx <- sample(1:nrow(df), size = 0.8 * nrow(df))
  train <- df[train_idx, ]
  test <- df[-train_idx, ]
  
  m_simple <- lm(Y ~ X1 + X2, data = train)
  m_ideal  <- lm(Y ~ X1 + I(X2^2), data = train)
  m_mixed  <- lm(Y ~ X1 + X2 + I(X1 * X2), data = train)
  
  test$pred_simple <- predict(m_simple, newdata = test)
  test$pred_ideal  <- predict(m_ideal, newdata = test)
  test$pred_mixed  <- predict(m_mixed, newdata = test)
  
  mse_simple <- mean((test$Y - test$pred_simple)^2)
  mse_ideal  <- mean((test$Y - test$pred_ideal)^2)
  mse_mixed  <- mean((test$Y - test$pred_mixed)^2)
  var_X1 <- var(test$X1)
  var_X2 <- var(test$X2)
  var_Y  <- var(test$Y)
  
  summary_table <- tibble(
    case = case_name,
    var_X1 = var_X1,
    var_X2 = var_X2,
    var_Y = var_Y,
    mse_simple = mse_simple,
    mse_ideal = mse_ideal,
    mse_mixed = mse_mixed
  )
  
  # Wykresy
  p1 <- ggplot(test, aes(x = X1, y = Y)) +
    geom_point(alpha = 0.5, color = "grey40") +
    geom_line(aes(y = pred_simple, color = "Simple"), linewidth = 0.8) +
    geom_line(aes(y = pred_ideal, color = "Ideal"), linewidth = 0.8) +
    #geom_line(aes(y = pred_mixed, color = "Mixed"), linewidth = 0.8) +
    labs(
      title = paste("Dopasowanie względem X1"),
      x = "X1", y = "Y", color = "Model"
    ) +
    theme_minimal()
  
  p2 <- ggplot(test, aes(x = X2, y = Y)) +
    geom_point(alpha = 0.5, color = "grey40") +
    geom_line(aes(y = pred_simple, color = "Simple"), linewidth = 0.8) +
    geom_line(aes(y = pred_ideal, color = "Ideal"), linewidth = 0.8) +
    #geom_line(aes(y = pred_mixed, color = "Mixed"), linewidth = 0.8) +
    labs(
      title = paste("Dopasowanie względem X2"),
      x = "X2", y = "Y", color = "Model"
    ) +
    theme_minimal()
  
  print(case_name)
  
  grid.arrange(p1, p2, ncol = 2)
  
  summary_table %>%
    kbl(caption = paste("Statystyki dla przypadku:", case_name)) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  full_width = FALSE)
}

plot_case_fit2("hidden_impact")
[1] "hidden_impact"

Statystyki dla przypadku: hidden_impact
case var_X1 var_X2 var_Y mse_simple mse_ideal mse_mixed
hidden_impact 9.894365 2.660213 1118.153 5.524909 0.0136601 0.0136601
Pokaż/Ukryj kod 
diagnose_model("hidden_impact")

Podsumowanie

Generowanie danych

Zdefiniowano “prawdziwą” funkcję generującą Y oraz globalny szum (ten sam poziom dla wszystkich scenariuszy), co umożliwia porównywanie MSE między przypadkami.

Rozkłady predyktorów: normalne (różne wariancje), jednostajne (wąskie/szerokie), wykładnicze, Poissona, t-Studenta, mieszane (w tym bimodalne) i kombinacje asymetryczne/ciężkoogonowe.

Modele i walidacja

Budowano modele dopasowane do struktury danych (proste vs „idealne” względem mechanizmu) i porównywano CV-MSE (K=5) wraz z przedziałami ufności.

Globalny szum utrzymywał wspólną skalę błędu, więc różnice w MSE wynikały głównie z rozkładów, postaci funkcjonalnej i obecności/siły Z.

Wpływ rozkładu na MSE

Rozkłady skośne i ciężkoogonowe (t-Student, wykładniczy, mieszane) zwiększają MSE i jego zmienność między foldami; modele są bardziej czułe na obserwacje skrajne.

Niedopasowanie postaci funkcjonalnej (np. pominięcie nieliniowości) podnosi MSE silniej niż sama zmiana rozkładu – to koszt uogólniania.

Standaryzacja i rozsądny wybór featurów stabilizują MSE, ale nie kompensują błędnej postaci modelu.

Zmienna ukryta

Z wprowadzono jako wspólną przyczynę, konstruowano proxy-Z (np. komponent główny).

Kluczowe obserwacje:

Specyfikacja modelu decyduje – poprawna postać potrafi „pośrednio” uchwycić wpływ Z i obniżyć MSE nawet bez jawnego Z.

Gdy wpływ Z jest dominujący, warto inwestować w dodatkowe zmienne/proxy/instrumenty; gdy dominuje błąd postaci, najpierw naprawiamy formę modelu.

Różnice MSE między modelem prostym a „idealnym” są czytelne i stabilne, co potwierdza powyższe tezy.

Epsilon progresu – realizacja

Diana Morzhak – pokazała, że wrażliwość MSE na strukturę danych jest użyteczną analogią do mechanizmów predictive coding

Mateusz Walo – przeprowadził causal approach: izolowanie efektu nieobserwowalnego czynnika, budowa proxy-Z, nadzorowanie i koordynowania prac i terminów realizacji projektu

Dominika Zydorczyk – usystematyzowała interpretację i komunikację wyników MSE dla różnych rozkładów/błędów, pokazując jak metryka odzwierciedla dopasowanie i kiedy myli modelarza.

Podsumowanie generalne

Zespół w krótkim czasie wykonał dużo pracy: od przemyślanego generatora danych, przez rzetelną walidację, po klarowne wnioski – każdy dowiózł swój epsilon progres, a razem osiągnęliśmy dojrzałe rozumienie, jak i dlaczego MSE zachowuje się tak, jak obserwowaliśmy.